function [criterion] = HPZ_MME_Analytical(observations,param,function_flag)

% number of observations
[num_obs,~]=size(observations);

% set prices from observation matrix
prices= observations(:,3:4);

% set observed bundle (x1,x2) from observation matrix
observed_bundle = observations(:,1:2);

% initialization of criterion matrix
criterion=zeros(num_obs,1);

if function_flag==1 % CRRA

    for i=1:num_obs

       % p = Max_Y / Max_X => p = p1 / p2
       p = prices(i, 1) / prices(i,2);
       
       min_x  = min (observed_bundle(i, 1), observed_bundle(i, 2));
       max_x  = max (observed_bundle(i, 1), observed_bundle(i, 2));
       rho  = param(2);
       beta = param(1);
       
       %% if Rho > 0 and Beat >=0
       if param(2) > 0 && param(1) >= 0
           
           %% p < 1/1+Beta
           if  p < 1/(1+param(1))
                   % x1* optimal:
                   x1_opt = ( ( (max_x^(1-rho)) + (1+beta)*(min_x^(1-rho)) ) / ( 1 + ( (1+beta)^(1/rho) ) * ( p^( (1-rho)/rho ) ) ) ) ^ (1/(1-rho));
                   % x2* optimal:
                   % x1*[(1+Beta)*p]^(1/Rho)
                   x2_opt = x1_opt * ( (1 + beta) * p ) ^ (1 / rho);
                   % compute the optimal expenditure, using optimal bundle
                   % (x1*, x2*)
                   criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
           end
           
           %% 1/1+Beta < p < 1+ Beta
           if p >= 1/(1+param(1)) && p <= (1+ param(1))
               % x1* optimal:
               x1_opt = ( ( (max_x^(1-rho)) + (1+beta)*(min_x^(1-rho)) ) / ( 2+beta ) ) ^ (1/(1-rho));
               % x2* optimal:x1*
               x2_opt = x1_opt;
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
           end
           
           %% p > 1 + Beta           
           if p > (1+ param(1))
               % x1* optimal:
               x1_opt = ( ( (max_x^(1-rho)) + (1+beta)*(min_x^(1-rho)) ) / ( (1+beta) + ( p / (1+beta) )^( (1-rho)/rho ) ) ) ^ (1/(1-rho));
               % x2* optimal:
               % x1*[p/(1+Beta)]^(1/Rho)
               x2_opt = x1_opt * ( p / (1 + beta ) )  ^ (1 / rho);
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
           end
                
       %% if Rho > 0 and -1 < Beat < 0 and x1_o >= x2_o:
       elseif param(2) > 0 && (param(1) > -1 && param(1) < 0)
           
           %% p < 1
           if  p < 1
               % x1* optimal:
               x1_opt = ( ( (max_x^(1-rho)) + (1+beta)*(min_x^(1-rho)) ) / ( 1 + ( (1+beta)^(1/rho) ) * ( p^( (1-rho)/rho ) ) ) ) ^ (1/(1-rho));
               % x2* optimal:
               % x1*[(1+Beta)*p]^(1/Rho)
               x2_opt = x1_opt * ( (1 + beta) * p ) ^ (1 / rho);
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
           end
           
           %% p = 1 
           if p == 1
               % x1* optimal:
               x1_opt = ( ( (max_x^(1-rho)) + (1+beta)*(min_x^(1-rho)) ) / ( 1 + ( (1+beta)^(1/rho) ) ) ) ^ (1/(1-rho));
               % x2* optimal:
               x2_opt = x1_opt * ( (1+beta)^(1/rho) );
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
           end
           
           %% p > 1          
           if p > 1
               % x1* optimal:
               x1_opt = ( ( (max_x^(1-rho)) + (1+beta)*(min_x^(1-rho)) ) / ( (1+beta) + ( p / (1+beta) )^( (1-rho)/rho ) ) ) ^ (1/(1-rho));
               % x2* optimal:
               % x1*[p / (1+Beta)]^(1/Rho)
               x2_opt = x1_opt * ( p / (1 + beta ) )  ^ (1 / rho);
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
           end
           
       %% if Rho == 0 and Beta >= 0
       elseif param(2)  == 0 && param(1) >= 0
           
           %% if  p < (1/(1+Beta))
           if  p <= 1/(1+ beta)
               x1_opt = max_x + (1+beta)*min_x;
               x2_opt = 0;
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
                       
           %% if p > (1+Beta)
           elseif p >= (1+ param(1))
               x1_opt = 0;
               x2_opt = max_x + (1+beta)*min_x;
               % compute the optimal expenditure, using optimal bundle
               % (x1*, x2*)
               criterion(i,1) = (prices(i,1) * x1_opt) + (prices(i,2) * x2_opt);
                         
           %% if  (1/(1+Beta)) <= p <= (1+Beta)
           elseif ( (p >= 1/(1+ param(1)) ) && ( p <= (1+ param(1)) ) ) && ( observed_bundle(i,1) >= observed_bundle(i,2) )
               criterion(i,1) = ( ( prices(i,1) + prices(i,2) ) * ( max_x + (1+beta)*min_x ) ) / (2+beta);
           end
           
       end
        
    end % end of loop over observations
else
    beta = param(1);
    A = param(2);
    for i = 1 : num_obs
        p = prices(i, 1) / prices(i,2);
        min_x  = min (observed_bundle(i, 1), observed_bundle(i, 2));
        max_x  = max (observed_bundle(i, 1), observed_bundle(i, 2));
        p_x = prices(i,1);
        p_y = prices(i,2);
        u_0 = -exp(-A*max_x) - (1+beta)*exp(-A*min_x);
        
        %% non-negative beta
        if beta >= 0 && A == 0 
            
            u_0 = max_x + (1+beta)*min_x;
            if p <= 1/(1+beta)
                criterion(i,1) = p_x*u_0;
            elseif p <= 1+beta && p >= 1/(1+beta)
                criterion(i,1) = ((p_x + p_y)/(2+beta))*u_0;
            elseif p >= 1+beta
                criterion(i,1) = p_y*u_0;
            end
            
        elseif beta >= 0 && -u_0 > 1+beta
            if p < - ((u_0 + (1+beta))/(1+beta)) %% Done
                criterion(i,1) = - (p_x * log(-(u_0 + (1+beta))))/A;
            elseif p >= - ((u_0 + (1+beta))/(1+beta)) && p < 1 / (1+beta)
                criterion(i,1) = p_x * ( log(-( (p+1)/(p*u_0) )) / A ) + ...
                                 p_y * ( log(- ( ((1+p)*(1+beta))/u_0 ) ) / A );
            elseif p >= 1/(1+beta) && p <= (1+beta) %% Done
                criterion(i,1) = (p_x + p_y) * (log(-(2+beta)/u_0) / A);
            elseif p > 1+beta && p <= - ((1+beta)/(u_0 + (1+beta))) %% Done
                criterion(i,1) = p_x * ( log(- ( ((1+p)*(1+beta))/(p*u_0) ) ) / A ) + ...
                                 p_y * ( log(- ( (p+1)/(u_0) )) / A );
            elseif p > - ((1+beta)/(u_0 + (1+beta))) %% Done
                criterion(i,1) = - (p_y * log(-(u_0 + (1+beta))))/A;
            end
            
        elseif beta >=0 && -u_0 <= 1+beta
            if p < 1/ (1+beta)
                criterion(i,1) = p_x * ( log(-( (p+1)/(p*u_0) )) / A ) + ...
                                 p_y * ( log(- ( ((1+p)*(1+beta))/u_0 ) ) / A );
            elseif p >= 1/(1+beta) && p <= 1+beta
                criterion(i,1) = (p_x + p_y) * (log(-(2+beta)/u_0) / A);
            elseif p > 1 + beta
                criterion(i,1) = p_x * ( log(- ( ((1+p)*(1+beta))/(p*u_0) ) ) / A ) + ...
                                 p_y * ( log(- ( (p+1)/(u_0) )) / A );
            end
            
        %% negative beta     
        elseif beta > -1 && beta < 0  && -u_0 > 1+beta %% Done
            if p < - ((u_0 + (1+beta))/(1+beta))
                criterion(i,1) = - (p_x * log(-(u_0 + (1+beta)) ))/A;
            elseif p >= - ((u_0 + (1+beta))/(1+beta)) && p <= 1 %% Done
                criterion(i,1) = p_x * ( log(-( (p+1)/(p*u_0) )) / A ) + ...
                                 p_y * ( log(- ( ((1+p)*(1+beta))/u_0 ) ) / A );
            elseif p > 1 && p <= - ((1+beta)/(u_0 + (1+beta))) %% Done
                criterion(i,1) = p_x * ( log(- ( ((1+p)*(1+beta))/(p*u_0) ) ) / A ) + ...
                                 p_y * ( log(- ( (p+1)/(u_0) )) / A );
            elseif p > - ((1+beta)/(u_0 + (1+beta))) %% Done
                criterion(i,1) = - (p_y * log(-(u_0 + (1+beta))))/A;
            end
            
        elseif beta > -1 && beta < 0  && -u_0 <= 1+beta %% Done
            if p <= 1
                criterion(i,1) = p_x * ( log(-( (p+1)/(p*u_0) )) / A ) + ...
                                 p_y * ( log(- ( ((1+p)*(1+beta))/u_0 ) ) / A );
            else
                criterion(i,1) = p_x * ( log(- ( ((1+p)*(1+beta))/(p*u_0) ) ) / A ) + ...
                                 p_y * ( log(- ( (p+1)/(u_0) )) / A );
            end

        %% beta equal to -1
        elseif beta == -1
            if p <= 1
                criterion(i,1) = p_x * max_x;
            else
                criterion(i,1) = p_y * max_x;
            end 
        end
        
    end %% end of for
    
end % end of CRRA and CARA

end % end of the MME function